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Abstract 

Calculation methods for turbulent duct flows are generalised for ducts with arbitrary cross- 
sections. The irregular physical geometry is transformed into a regular one in computational 
space, and the flow equations are solved with a finite- volume numerical procedure. The turbulent 
stresses are calculated with an algebraic stress model derived by simplifying model transport 
equations for the individual Reynolds stresses. Two variants of such a model are considered 
in the present study. These procedures enable the prediction of both the turbulence-driven 
secondary flow and the anisotropy of the Reynolds stresses, in contrast to some of the earlier 
calculation methods. Model predictions are compared to experimental data for developed flow in 
triangular duct, trapezoidal duct and a rod-bundle geometry. The correct trends are predicted, 
and the quantitative agreement is mostly fair. The simpler variant of the algebraic stress model 
procured better agreement with the measured data. 


*On leave from University of Lagos, Lagos, Nigeria. Work funded under Space Act Agreement C99066G. 
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NOMENCLATURE 


B 

half-width of duct 

Cl 

constant in algebraic stress model 

Cn,C e \,C(2 

k — e turbulence model constants 

D 

determinant of the Jacobian of the transformation 

f 

near-wall proximity function in algebraic stress model 

G 

rate of production of turbulent kinetic energy 

H 

half-height of duct 

J 

Jacobian of the transformation y : = y J (x J ) 

k 

turbulent kinetic energy 

L 

mixing length 

P 

pressure 

9* 

stress or flux of variable $ 

s 

distance from point in flow to a point on wall 

s+ 

source term for dependent variable $ 

u b 

bulk velocity 

Uo 

maximum velocity 

u. 

friction velocity 

Ui or U,V,W 

cartesian velocity components; streamwise, lateral, vertical 

h? 

hr 

hr 

Reynolds normal stresses in cartesian directions 

UiU2,UlU3,U 2 U 3 

Reynolds shear stresses (2 primary and 1 secondary) 


curvilinear coordinates 

y j 

cartesian coordinates 

< y 2 > 

square of the average distance from a point to all walls 

a, P, 7 

constants in algebraic stress model 

«*• 

metric coefficients of transformation 

Sii 

Kronecker delta 

€ 

rate of dissipation of turbulent kinetic energy 

6 

angle 

K 

von Karman constant 

pt 

turbulent eddy viscosity 

p 

density 

<7<f, 

turbulent Prandtl/Schmidt number for 

T 

wall shear stress 

$ 

General representation of dependent variable 

Superscripts 

i 

longitudinal direction 

2,3 

vertical and lateral directions 


2 



1 INTRODUCTION 


Secondary motion of PrandtPs first kind is encountered in curved passages, where centrifugal forces 
act at an angle to the streamwise direction. It is pressure-driven, with magnitudes ranging from 
10-40 % of the bulk streamwise velocity. It usually leads to a first-order perturbation of the main 
flow properties, and there is no real difficulty in calculating it. On the other hand, secondary 
motion of PrandtPs second kind is turbulence-driven and is believed to result from the anisotropy 
of the turbulent stresses. It is associated with internal or external turbulent flows in the presence of 
streamwise corners. Its magnitude has been observed to be mostly between 1-4 % of the bulk mean 
velocity in most straight ducts with non-circular cross-section, but the effects on wall shear stress 
distribution and heat transfer rates are quite significant. Some form of Reynolds stress modelling 
is usually required in order to calculate this type of secondary motion. 

Probably the first qualitative data on turbulence-driven secondary motion in non-circular ducts 
were presented by Nikuradse[l], using flow visualisation techniques. Actual quantitative measure- 
ments were presented by Hoagland[2] aboout 30 years later. Several detailed measurements have 
subsequently been presented. A comprehensive review of experimental data for ducts with regular 
cross-section can be found in Demuren and Rodi [3], Additional data for irregularly-shaped chan- 
nels are reported by Rodet [4], Aly et al [5], and Seale [6], for trapezoidal, triangular and simlated 
rod-bundle-type cross-sections, respectively. 

The first calculation of secondary flow in straight non-circular ducts was due to Launder and 
Ying [7] (hereafter denoted LY). Their method is based on approximating the Reynolds stresses with 
algebraic expressions by simplifying the corresponding transport equations proposed by Hanjalic 
and Launder [8] This model has subsequently been employed, sometimes with modifications in 
numerous studies [9-12] to calculate secondary flows in ducts with regular as well as gemetrical cross- 
sections. In all these studies, it has been possible to predict the magnitudes of the secondary motion 
fairly well, but always at the expense of the Reynolds stress distributions. The anisotropy of the 
secondary normal stresses are usually underpredicted by an order of magnitude. This discrepancy 
has been discussed by various authors, including Demuren and Rodi [3], who showed that certain 
terms which are neglected in LY model could in fact be larger in magnitude than terms which 
are included. In order to avoid this shortcoming, Seale [6] in the calculation of developed flow in 
a rod-bundle-type configuration, utilises an algebraic model for the streamwise vorticity source, 
based on the correlation of measured turbulence data. Very good agreement between measured 
and computed velocity and wall shear stress distributions are reported. The prospects of such 
a method as a general simulation technique for arbitrary geometries is uncertain. The reported 
calculations are perfomed on a fairly coarse cartesian grid, using a method which approximates the 
curved wall with steps. The influence of such a boundary treatment on the total flow is uncertain, 
but it doubtful that it can produce a true picture of the flow near the curved wall. A truly general 
computational method should employ a boundary-fitted coordinate system. 

Naot and Rodi [13] (hereafter denoted NR) calculated flows in rectangular ducts using an 
algebraic stress model developed from the Reynolds stress model of Launder Reece and Rodi [14]. 
The method was subsequently extended by Demuren and Rodi [3] (hereafter denoted DR). To 
put the various models in perspective. LY model is derived from the Reynolds stress model of 
Hanjalic and Launder [8]. The latter is a precursor to the more complete Reynolds stress model 
of Launder, Reece and Rodi [14]. The LY model neglects all terms involving secondary velocity 
gradients. NR model is based on the more complete Reynolds stress model and approximates terms 
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involving secondary velocity gradients by eddy viscosity type relations. DR model is also based on 
the same Reynolds stress model, but uses the full derivation. In all the models the algebraic stress 
relations are derived from the partial-differential equations governing the Reynolds stress transport 
by neglecting the convection and diffusion terms or relating them to those of the turbulent kinetic 
energy. This approximation is adequate in straight wall bounded flows, but would be inadequate 
for free shear flows or flows with high mean shear rates in which the convection and diffusion terms 
assume more significant roles. As may be expected, DR model is more difficult to implement than 
the other two. It is also much more sensitive to the solution procedure, because of the much stronger 
coupling between the Reynolds stresses. Both NR and DR model produce correct simulation of 
mean flow and Reynolds stresses in developing and fully developed flow in square and rectangular 
ducts [3,15]. In fully developed flow situations the secondary flow level predicted with the former 
agree better with experimental data. 

In the present study, NR and DR models are generalised for any curvilinear coordinate system 
and then applied to predict developed turbulent flow in straight ducts with arbitrary cross-sections. 
The governing equations are transformed into a general curvilinear coordinate system, but the 
cartesian velocity components are retained as dependent variables, following the practice of Rhie 
and Chow [16]. 

2 MATHEMATICAL MODEL 

The time-averaged, three-dimensional, steady-state equations governing turbulent flows may be 
expressed in curvilinear tensor, for a general property $ as : 

-§2 \D(pUJt - »*.)<£) = DSt, (1) 

where x 3 is a general curvilinear coordinate system related to the cartesian coordinate system 
y J by the transformation y J = y J (x J ). 

D is the determinant of the Jacobian of the transformation; i.e., D = | J|, with 



tty, also called the metrics, are the elements of the Jacobian J 1 of the inverse transformation, 
x j = x j (yi), i.e., oc) = 

U m represents the cartesian velocity components. 

q^m and are respectively fluxes/stresses and source terms for the corresponding $. They 
may be expressed as : 
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Conservation Equation 

$ 


s* 

Continuity 

1 

0 

0 

Momentum Components 

Ui 

Pf>mi PUm^i + fi + fi x T a i) 

0 

Turbulent kinetic energy 

k 

<*k clx nQ!m 

G - pe 

Dissipation rate of k 

€ 




where, is the Kronecker delta 

and, G is the rate of production of turbulent kinetic energy by the interaction of turbulent 
stresses with mean shear, and is given by : 


g = MiK + f*W) 

The turbulent eddy viscosity fit is given by the standard k — e model as 

k 2 

Ht = c^p^- 

The empirical constants are: 


( 2 ) 

( 3 ) 


CfX — 1.44 , c e2 = 1.92 , a k = a t — 1.225 (4) 

In the above, we have used the Einstein summation convention for repeated indices to write in 
compact form. Thus, for the 3-dimensional flow under consideration, 

°ti>? = <4&i + 4&; + ( 5 ) 

The equations can be simplified if we consider only straight channels with constant cross-section 
in the streamwise direction (y 1 = x 1 ). Then, we can use parabolic assumptions to decouple the 
pressure gradient in the streamwise direction from those in the cross-plane, and also neglect the 
derivatives with respect to x 1 or y 1 , in all but the convection terms. They will of course be zero 
throughout in fully developed flow, and the equations become 2-dimensional. 

The major difficulty in solving the set of equations is to determine the distribution of the 
Reynolds stresses —puiuj. Following Demuren and Rodi [3,15] , we can write the algebraic equation 
for the Reynolds stresses as : 


~(i - a )i u i u iffii& a r + u j u i§^& a r] 

+ /? 

- [§£&<*? + !pr<*r] - £ [ciU^r- + |(1 - a - (3 - c x ) Sijk] = 0 (6) 

Equation (6) shows very strong coupling between the Reynolds stresses, and in order to stabilise 
the solution process only the secondary stresses u\ , u\ and u 2 us are obtained from it. The primary 
normal and shear stresses are calculated from the turbulent kinetic energy and standard eddy 
viscosity relations. The resulting set of algebraic stress relations are:- 


5 







( 7 ) 

( 8 ) 
( 9 ) 


u\ = 2k — u\ — u\ 

-to = ^f|W 

-TO = 

u 2 = + P + C 1 — 1 ) + P^l ^2§^ a T 

~ U2“3{(1 - a )fi^m a T “ 

“ 7*f|W] / [1 + cfH 1 ~ a ~ P)§^ a T] ( 10 ) 

u 3 = Z\e\{° t + P + c i ~ 1) + P^H^^k^T 
- U 2 U 3 {(1 - a)^~^a^ - P^ik&T} 

~ Tf*JpWl / [1 + ^(1 - a - $)§gkc?) (11) 

“2®3 = c7e[^( UlU 3^p^ a r + ^1^2 §^ a T) 

- {(1 - <x)u\ - Pul + Tf fc}|p:ar - {(! - a )«l - 0«1 + ^fpW] (12) 

In equations (10 - 12) wall proximity effects are accounted for by making the empirical coeffi- 
cients functions of a dimensional distance from the walls, thus : 


a = 0.7636 - 0.06 /; fi = 0.1091 + 0.06 /; 7 = 0.182; ci = 1.5 - 0.50 / (13) 

where f is the wall-proximity function given by : 


/ = 


n . 

< y* > ’ 


<y 2 > = 


2 T 2 * 

w Jo 1? 


(14) 


C 13/2 

L = is the length scale of the turbulent motion at the grid node, and < y 2 > 1 / 2 is 

the average distance from all walls. 

In the NR model, equations (10 - 12) are further simplified by setting the denominator to unity 
and approximating terms involving secondary velocity gradients by eddy-viscosity-type relations 
as: 


£[TO{( 1 - «)J|W - /*f|W} - ( 15 ) 

^[u^uE{(l - a)Jr^a™ - P^£.a™} - lk^£a™] ~ ( 16 ) 

cTe [{( x “ a ) u 2 ~ P u 1 + 
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- {(1 - a)u | - /3u f + 7 + §^ a T) ( 1 7 ) 

This practice has a stabilising effect on the solution procedure, since they can now be treated 
implicitly as diffusion terms in the momentum equations for the secondary velocity components. 
An analysis of the implications of these approximations, performed by Demuren and Rodi [3] for 
a square duct flow, has shown that the major effect is on the secondary shear stress calculation, 
of which equation (17) is a dominant part. That term is essentially reduced by a factor of 2 - 3. 
However, the secondary shear stress acts a sink of the streamwise vorticity, so that a reduction 
in its magnitude would lead to an increase in the secondary motion. This explains why the NR 
model usually produces stronger secondary motion than the DR model [3]. The terms in equations 
(15) and (16) are small compared to the rest of the terms in equations (10) and (11), so that the 
effects of approximating them is minimal. If the calculations are performed with the secondary 
velocity terms in equations (15) - (17) neglected completely, as is done in LY model and derivatives 
thereof, the secondary motion will grow indefinitely and no stable solution can be obtained. In this 
case, a stable solution can only be obtained through an ad-hoc reduction in the anisotropy of the 
secondary normal stresses, which is the source of streamwise vorticity. 

3 SOLUTION PROCEDURE AND BOUNDARY CONDITIONS 

The differential equations to be solved are parabolic in the streamwise direction so that influences 
only travel from upstream to downstream. An efficient once-through forward marching procedure 
can thus be used, starting from given inlet conditions. The procedure is essentially an incorpora- 
tion of a non-staggered grid arrangement into the SIMPLE algorithm of Patankar and Spalding 
[17]. Where the grid is not orthogonal, cross-derivative terms lead to the formation of a 9-point 
computational molecule for the 2-dimensional plane. In the present study, these terms are treated 
explicitly in the equations for the velocity components. This practice proved not to be so robust 
in the solution of the pressure correction equation, which is the most critical for stability and 
convergence. We therefore employ a 9-point pressure-correction equation which is an extension of 
Stone’s strongly implicit scheme as proposed by Peric [18]. A comparative study [19] of various 
solution schemes, including coupled and decoupled methods, for two-dimensional flow problems on 
a non-orthogonal grid shows the present arrangement to be the most effective. 

Four types of boundaries are usually encountered, at which conditions must be specified for 
the dependent variables. At inlet planes, known (or simulated) conditions are specified for all 
variables; usually uniform values for the primary velocity and the turbulent quantities, and zero for 
the secondary velocity components. No specific conditions are required at outflow boundaries, since 
the equations are parabolic. Along a symmetry plane, the velocity component normal to the surface 
is set to zero. Further, the gradient of all other variables is zero in the direction normal to the 
surface. At solid walls, the velocity component normal to the surface is also set to zero, but rather 
than carry out the solution all the way down to the wall, we employ the wall-function method [3] to 
specify the values of the velocity components parallel to the surface, and the turbulence quantities, 
along the first line of grid nodes away from the wall. These grid nodes are usually arranged to be in 
the fully turbulent region, where the logarithmic law of the wall may be assumed to be applicable. 
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4 RESULTS AND DISCUSSION 


4.1 Test Cases 

Three test cases are considered in the present study. These coincide with flow in ducts with 
cross-sections requiring non-rectangular grids, for which experimental data could be found in the 
literature. In each case, the channel cross-section is constant in the streamwise direction, and the 
flow is nominally fully-developed. 

Case I is the flow in an equilateral triangular duct (fig. 1) studied experimentally by Aly 
et al [5] at Reynolds number of 53,000, based on the bulk velocity and the hydraulic diameter. 
Measurements were made in one-sixth of the duct cross-section at a plane 133 diameters from the 
inlet. Case II is the flow, at a Reynolds number of 240,000, in a trapezoidal channel with two sides 
parallel and two slanting sides at angles 75° and 105° (fig. 2). Rodet [4] presents measurements 
of the velocity distribution and the Reynolds stresses at a cross-stream plane 50 diameters from 
the inlet. Judging by the review of duct flow data, this would barely be long enough for full flow 
development. The cross-section of the channel for case III is shown in fig. 3. It was studied 
experimentally by Seale [6] to simulate approximately flows encountered in rod-bundle geometries. 
Measurements were made in a symmetric quadrant at a Reynolds number of 200,000 and 120 
diameters. 

4.2 Computational Details 

In the calculations, advantage is taken of the inherent symmetry of the flow to reduce the size 
of the computational domain. Thus, we compute only a third of the triangular duct half of the 
trapezxoidal duct, and a quarter of the rod- bundle-type geometry. In each case, the grid is generated 
algebraically using transfinite mapping of the boundary points. Where necessary, a third-order 
smoother based on Newton’s divided difference formula is applied. Although orthogonality is 
desirable, it is not a requirement of the present solution procedure. The sole constraint applied is 
that the normal distance of the first set of grid nodes, nearest to the wall should lie in the range 
y+ - 30 - 200. 

Since all the present test cases apply to fully-developed flow, the grid distribution in the stream- 
wise direction is unimportant. The practice adopted here is to take an infinitely large step size, and 
then use a relaxation method to obtain a converged solution. An underrelaxation method with a 
maximum value of 0.5 is applied to all variables. A typical computation on a 40 x 40 grid required 
about 600 iterations to converge to a normalised residual norm of 10~ 6 when the NR model is used. 
DR model calculations require 200-300 iterations more. Grid dependency tests are performed by 
doubling the total number of grid nodes, i.e., an increase by a factor of \/2 in each direction, and 
comparing the results. These tests show that a 26 x 26 grid is sufficient for case I. Case II requires 
a 18 x 36 grid, and case III a 40 x 40 grid. 

4.3 Comparison with Experimental Data 
4.3.1 Triangular Duct 

The measured wall shear stress distribution is compared to present predictions with the k — € , NR 
and DR models in fig. 4. Both the NR and DR models capture the trend of the data towards more 
uniformity, though the former has the better agreement. The k — e model predicts much wider 
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variation in wall shear stress. These can be explained with the results of the predicted secondary 
flow presented in figs. 5 and 6. This flow moves towards the corners along the corner bisector 
and returns along the walls towards the wall bisectors. Whilst the NR and DR models predict the 
correct trend, the k — e model produces no secondary motion. This motion leads to a redistribution 
of the wall shear stress by moving high velocity fluid from the central part of the duct to the 
corners thereby increasing the shear. In contrast, the shear stress near the wall bisector is reduced 
by the slower moving fluid convected along the wall. The influence of this secondary motion on the 
streamwise velocity distribution is shown in fig. 7. The experiments show strong bulging of the 
velocity contours towards the corners. The predictions mirror this effect to varying degrees. The 
k — e model predicts nearly circular contours (not shown). Fig. 8 shows similar comparisons for 
the turbulent kinetic energy. The normal stresses are compared in fig. 9. Both the NR and DR 
models predict the correct degree of anisotropy. 

4.3.2 Trapezoidal Duct 

The present predictions with both NR and DR models show that there are two eddies in each 
corner. The one stradling the longer side is much larger, as may be expected. Fig. 10 compares 
the predicted contours of the streamwise velocity component to the measurements of Rodet [4]. 
NR model predicts the bulging of the contours towards the corners to the correct degree, though 
there is some deviation between the actual levels along the centreline. DR model predictions agree 
better with the magnitudes along the centreline, but produces too little bulging of the contours 
towards the corners, due to its characteristic of producing weaker secondary motion. Predictions 
of the turbulent kinetic energy, k, are compared in fig. 11. The effect of the secondary motion 
in distorting the distribution of k is much stronger. Correspondingly, the differences between the 
predictions are magnified. Again NR model shows the correct trends overall, but DR model has 
better agreement with centreline values. 

The source of streamwise vorticity is the anisotropy of the secondary normal turbulent stresses. 
A measure of this anisotropy is shown in fig. 12. The present prediction (NR model) compares 
quite well with the data. The predictions by Nakayama et al [12], which is characteristic of those 
based on the LY model show magnitudes which are at least an order of magnitude too low. This 
observation has, of course, been widely discussed in the literature without isolating the true cause. 

4.3.3 Rod-Bundle Geometry 

Case III represents the most complicated geometry of all the test cases computed in this study. 
There are 5 corners, all different from each other. We have combinations of intersections of straight 
walls, curved walls and symmetry planes. The geometry certainly poses a severe test on the pre- 
dictive ability of the present mathematical procedure and turbulence model. There are, however 
significant differences from the usual rod-bundle geometries, which contain no wall/wall intersec- 
tions. Figure 13 compares the wall shear stress distributions predicted with the three turbulence 
models to measured data of Seale [6]. The data shows nearly uniform wall shear stresses, away 
from the corners. This effect is, of course, due to the equalizing action of the the secondary motion. 
The k - e model cannot simulate this effect since it does not produce any secondary motion. DR 
model does produce a more uniform distribution, but the secondary motion still appears to be 
too weak. NR model predicts the correct level of shear stress, but this has a wavy distribution, 
not apparent in the measurements. It is not certain if the waviness is physical, but the results 
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have been reproduced with the NR model, using a different numerical procedure, so it does not 
appear to be numerical in origin. The measurements do not appear to be sufficiently detailed to 
provide definitive answer. Examination of the lateral velocity data, presented in fig. 14, shows 
very strong periodic variations along the verticals. The predictions do simulate this behaviour. As 
usual, the NR model performs better. The peaks in the predicted wall shear stress profiles are 
seen to coincide with the maxima in the velocity component directed towards the vertical wall. 
Thus the data suggests the action of the secondary motion near the vertical wall could produce 
miltiple peaks in the shear stress distribution as predicted. A similar effect would be expected on 
the turbulent kinetic energy. These measurements represent the most detailed set presented by 
Seale [6], and the good agreement with the predictions provides the confidence that the correct 
behaviour is being simulated. Fig. 15 compares measured and predicted vertical velocity profiles. 
Again we see a periodic variation in the predictions, but the measured data points are too scanty 
to discern any trend. Considering the measured data, Seale claims large continuity errors ( 0.19 
[/*) along horizontal lines. Close to the horizontal wall, the imbalance was as high as 1.78 f/*, 
which is considerably larger than the maximum secondary velocity measured. The presented data 
certainly do not provide an accurate representation of the flow, and have been included here only 
for completeness. Since the lateral velocity field is well predicted, continuity would ensure that the 
vertical velocity field is too. Figure 16 shows the velocity vectors and the stream function contours 
predicted with the NR model. The eddy formation appears to be quite complex. There are at least 
two eddies at each wall/wall juction, and the wall/symmetry plane junctions each contain one eddy. 
The two large central eddies appear to be in the process of breaking up, as has been observed in 
channels with large aspect ratios. The predictions with the DR model are similar, but with reduced 
strengths. 

Predicted contours of the streamwise velocity component and the turbulent kinetic energy are 
compared to measured data in figures 17 and 18, respectively. The distortion of the contours due to 
the secondary motion is well predicted, but the quantitative agreement is only fair. The predicted 
components of the turbulent intensities along the vertical centreline are compared to correlations 
of the measured data in fig. 19. The agreement is fairly good, considering that there could be 
uncertainties of about ± 29 % in the correlation coefficients [6], 

It is pertinent to mention that the wall-function method used here precludes adequate resolution 
of the secondary velocity profiles very close to the wall. However, there is as yet no consistent 
Reynolds stress model for this near- wall region. Mansour et al [20], have discussed the inadequacy 
of the standard Launder, Reece and Rodi [14] model in this region. Rodi [21] presents a review of 
current attempts to bridge this gap. Preliminary calculations by the present author for a developed 
square duct flow, using such a near-wall model, shows that the computed mean flow and turbulence 
quantities are essentially the same as those obtained with the present model in the region more 
than 1 % of the duct width, away from the wall. However, in consonance with other attempts [21], 
the predicted components of the Reynolds stresses do not agree with the limiting two-dimensional 
behaviour of turbulence in the region very close to the wall. It appears that major changes to the 
pressure-strain modelling would be required to overcome this problem. 

5 CONCLUDING REMARKS 

Calculation of turbulent flow in complex geometries have been presented. The performance of 
three turbulence models are compared to measured data for three different channel configurations. 
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The standard k — e turbulence model does not predict any secondary motion of Prandtl’s second 
kind, and thus fails to predict important flow features. Although the NR model can be seen 
as an approximation of the DR model, it procures better quantitative agreement with measured 
secondary flow data. It thus enables the correct degree of redistribution of the primary velocity 
and the turbulence fields that takes place in complex channels to be predicted. Computations with 
the DR model also show the correct trends, but the magnitudes of the secondary flow is mostly 
underpredicted. Previous analysis had shown that differences in the treatment of the secondary 
shear stress are responsible for the different strengths of secondary motion predicted. Unfortunately, 
this stress component is the most difficult to measure, and it has received the least attention in 
turbulence model studies. More attention will have to be given to this component if there is to 
be significant progress in turbulence modelling of flow in channels with complex configuration. 
Coupled with this is a need for better near-wall modelling. The widely used wall-function approach 
now appears inadequate, for accurate simulation of flow in such complex geometries. Work is 
progressing on such a model. It must ensure that the various components of the Reynolds stresses 
exhibit the correct limiting behaviour as the wall is approached. Since detailed experimental data 
is limited, we have to make use of data generated by direct simulation of channel flows, such as 
those by Kim et al [22], and analysed by Mansour et al [20] in the light of popular Reynolds stress 
modelling concepts. 
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FIGURE 2. - FLOW CONFIGURATION AND GRID 
DISTRIBUTION FOR TRAPEZOIDAL DUCT. 
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FIGURE 3. - FLOW CONFIGURATION AND GRID DISTRIBUTION FOR 
SIMULATED ROD-BUNDLE GEOMETRY. 
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FIGURE 9. - TURBULENT INTENSITIES ALONG CENTRELINE. 
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